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Abstract 

An adaptive implicit finite difference method with non-uniform timesteps for 
solving the fractional diffusion equation in the Caputo form is proposed. The 
method allows one to adapt the size of the timesteps to the behaviour of the 
solution in order to keep the numerical errors small without the penalty of a huge 
computational cost. The method is unconditionally stable and convergent. In 
fact, it is shown that consistency and stability implies convergence for a rather 
general class of fractional finite difference methods to which the present method 
belongs. The computational advantages of the present adaptive method against 
fixed step methods are illustrated by solving the problem of the dispersion of a 
flux of subdiffusive particles stemming from a point source. 
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1. Introduction 

Fractional calculus, a very old field of Mathematics dating back to the time 
of Leibnitz has recently become a research area of growing interest due 
mainly to its surprisingly broad range of applications in Physics, Engineering, 
Chemistry, Biology, Economics [1, U 0, @, @ 0, B @1 ■ In particular, fractional 
calculus is a key tool in the study of some anomalous diffusion processes, which, 
for example, are abundant in biological environments where the presence of 
traps and obstacles often leads to mean square displacements that grow sub- 
linearly with time 0, [Tol 11|. A well-known model for describing this kind of 



subdiffusion processes is the so-called Continuous Time Random Walk model in 
which waiting time distributions between successive steps have a power-law tail. 
A convenient property of this mesoscopic model is that it leads to macroscopic 
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fractional diffusion equations that describe how the concentration of walkers 
evolves in space and time [3, [HI, [l3j] . 

Some analytical methods of solution of these equations are known (method 
of ima ges, separation of variables, integral transform methods,. . . ) 0, [a, 14, 15 



3 JJZ , 18 1 . However, as is also the case for non- fractional problems, often it is 



not possible to find an (at least useful) analytic solution and one has to resort 
to numerical methods. The proposal, development, and analysis of numerical 
methods to solve fractional differential equations is at present a quite active 
field of research, and many and varied methods have been considered Il9l . I20L 

HL [13, ES [S3, [11, m, [13, [il, [2^, [13, HH [S3, m, [H, m, [13, HE [1M S3, SiL Elf. 

Among these methods, finite difference methods are particularly convenient, 
and consequently have been extensively studied. However, in almost all the 
cases considered the timesteps are constant (a recent exception is the work by 
Polubny et al. [3!|). Certainly, in this way these methods are simpler. On the 
other hand, methods that allow variable timesteps have the great advantage that 
they can be converted into adaptive methods in which the size of the timesteps 
is chosen according to the behaviour of the solution. For example, one can 
choose small timesteps when the solution changes quickly in order to avoid the 
typical problem of methods with constant timesteps of passing over, or scarcely 
sampling, those regions where the solution has this behaviour. Besides being 
more reliable, adaptive methods are usually faster because they can employ 
large timesteps whenever the solution changes smoothly. In particular, one 
could dynamically adjust the size of the timestep so that the error is smaller 
than a prefixed value. The aim of this paper is to extend this procedure to 
fractional diffusion equations. 

The use of adaptive finite difference methods to fractional equations is espe- 
cially relevant because, for these equations, the number of operations required 
to calculate the numerical solution at the nth timestep scales as n 2 . In compar- 
ison, this number simply scales as n for normal equations. The n 2 behaviour 
is due to the requirement of using the values of the numerical solution for all 
the previous times at which the solution was calculated (because the fractional 
derivatives are non-local integro-diffcrential operators), which means that the 
number of calculations is roughly proportional to $2 i=1 i ~ nr. This makes 
these methods slow and hugely memory demanding when n is large. The prob- 
lem is so acute that a great deal of effort has been devoted to overcoming it. 
Two main approaches have been explored. The most obvious one consists of in- 
creasing the ability of the numerical method to evaluate the non-local operators 
with large timesteps without loss of accuracy 43|, |44[ . A less straightforward 
approach is that of the nested meshes method which is based on some scaling 
properties of the fractional derivatives (45[ . The number of operations in this 
method scales only as nlnn. It should be noted that these two approaches are 
compatible with the adaptive method with non-uniform temporal meshes to be 
discussed in this paper. 

The stability of a finite difference method is, from a practical point of view, 
the key property one has to assess because it is not useful numerically if one does 
not know under what circumstances the method is unstable. The stability of the 
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present adaptive method is analyzed by means of the von Neumann (or Fourier) 
procedure particularized for fractional diffusion equations [21 , 22|, 23|, 0, [3l[ . 



It is found that the method is unconditionally stable regardless of the temporal 
non-uniform mesh employed. Besides, it is proved in a rather general form that 
under relatively weak conditions the consistency and stability of a fractional 
finite difference method implies its convergence, i.e., it implies that the solution 
of the continuous equation is recovered from the finite difference solution when 
the size of the discretization mesh goes to zero (46. 47. 



We shall consider the fractional diffusion equation in the Caputo form 

du = F (la) 



with 
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d~> 1 /"* 1 dy(r 



and where 

^^w^)L dT w^^ 0<7<1 (2) 

is the Caputo fractional derivative and F(x,t) is a given source term [l9|. It 
should be noted that our procedure can be extended to other equations with 
other terms (fractional Fokker-Planck equations, fractional diffusion-wave equa- 
tions where 1 < 7 < 2) or even with other fractional operators such as the 
Riemman-Liouville derivative. In some diffusion problems u(x, t) represents 
the probability density of finding a particle at x at time t. When the particle 
starts at zero at time zero, the solution u(x, t) of du = is the propagator 
(or Green's function) with u(x, 0) given by Dirac's delta function, S(x), and 
(x 2 ) ~ 2Kt y /T(l + 7) is the corresponding mean square displacement of the 
particle for large t. For this reason the parameter K is sometimes called the 
anomalous diffusion coefficient and 7 the anomalous diffusion exponent. The 
operator d 1 j dt 1 is the usual first order derivative when 7 = 1 

The fractional diffusion equation is also usually written in terms of the frac- 
tional Riemman-Liouville derivative [!, HBj]. From a practical perspective, the 
two representations are different ways of writing the same equation as they are 
equivalent under fairly general conditions [25| . However, the numerical methods 
obtained for each representation are formally different, which leads to different 
ease of use and efficiency of the numerical algorithm. For example, in [25| it 
was proved that the method of Gorenflo et al. [2(| for the diffusion equation in 
the Caputo form and the method of Yuste and Acedo in [23| for the diffusion 
equation in the Ricmann-Liouvillc form were equivalent only if the Grunwald- 
Letnikov discretization (or BDF1 formula) is used in both methods. 

The difference scheme used in this paper is obtained by discrctizing (i) the 
Laplacian with the three-point centred formula and (ii) the Caputo derivative 
with a generalization of the LI formula to non-uniform meshes. It is second- 
order accurate in the spatial mesh size and first-order accurate in the timestep 
size. For a uniform temporal mesh, it becomes the numerical scheme discussed 
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by Liu et al. [30( and Murio 31j. We choose this scheme because it is suitable 
for being transformed into an adaptive method with non-uniform timesteps and, 
also, because it allows an easy consideration of discontinuous solutions (see 
Section [5]) • For example, this is not true for schemes based in the Griinwald- 
Letnikov discretization. 

The paper is organized as follows. In Section [5] the finite difference scheme 
with non-uniform timesteps is obtained. In Scction[3]thc stability of the method 
is analyzed by means of the von-Neumann stability method. Its convergence is 
proved in Section 2) In Section [5] the method is applied to the problem of 
evaluating the dispersion of a constant flux of subdiffusive particles stemming 
from a given point source. The paper ends with some conclusions and remarks 
in Section [6] 



2. Fractional difference algorithm with non-uniform timesteps 

As usual in finite difference methods, one starts by considering a mesh in 
the space-time region where one wants to obtain the numerical estimate C/j 1 "^ 



of the exact solution u(xj,t m ) = uj m; , (xj,t m ) being the coordinates of the 
(j,m) node of the mesh. Next, one replaces the continuous operator d of the 
equation du = F one has to solve by a difference operator 5 and a truncation 
error R(x, t): du = du + R. For example, in this paper we replace the continuous 
operators that define d in (flbf by 

^u{x,t) = 5?u{x,t) + Rt{x), (3) 
d 2 

-^u(x,t)=S 2 x u(x,t) + R x {t), (4) 

where 8] and S 2 are the corresponding difference operators. This way one gets 

[8] - KSl] u(x, t) = F(x, t) + R(x, t) (5) 

where R{x 1 t) = KR x (t) — Rt(x). In this case 5 = 6] — K8 2 . Neglecting the 
truncation term, one gets a difference equation SU = F whose solution leads to 
the finite difference estimate of the exact solution u(x, t) at the mesh points. For 
a given operator d one can consider many different difference operators 5, and 
hence many different finite difference methods for solving the finite difference 
equation du = F 0, 0, S3 . 



Wc here assume that the spatial size of the mesh Xj + \ — Xj = Ax is con- 
stant, and discretize the Laplacian operator d 2 jdx 2 by means of the three-point 
centred formula so that 

2 u(x ]+ i,t) - 2u(xj,t) + u(x j+1 ,t) 

Wji f J- (Ax) 2 

with Rt(xj) = O(Ax) 2 . For the fractional Caputo derivative we choose a dis- 
cretized operator 6] that is a generalization of the LI formula [H for non-uniform 
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meshes (see the Appendix): 

^ n— 1 

5?u(x,t n ) = _ ^ T ™1 [ufatm+l) - u ( x ^m)] > ( 7 ) 
m=0 

Rt{n) being of order i* -7 maxo< m <„-i (t m +i — i m )- Therefore the truncation 
error is 



R{Xj,t n ) = 0(Ax) 2 +t 1 -'tO 



max (i m+ i - t m ) 

0<m<n-l 



(8) 



which goes to zero when the spacing of the time-spatial mesh goes to zero. This 
means that the present method is consistent [471 ]. 

Therefore, introducing and ([7]) into S = 5] — KS 2 , the finite difference 
equation SU = Fwe have to solve becomes, after multiplying it by t n — i„_i, 
dU = F, i.e., 



m=0 



where 





= r(2- 


7j (Ax) 2 ' 


7> (7) 

m,n 


= {t n ~ 


+ A7t;(7) 

I'n—l) J - m ,ni 


T (7) 

0,1 


= (*!" 




7^(7) 
m,n 


(*n- 










= (*«- 





m < n — 1 , 



(9) 

(10) 

(11) 
(12) 

(13) 

(14) 



From here on we use the convention that a tilde over a symbol stands for that 
symbol multiplied by t n — i n _i. Reordering ([9]) and taking into account that 
T'^ 7 _''i n = 1 [see (|13l) ] we get the finite difference scheme we were looking for: 

- S n U? +1 + (1 + 25„)U? - S n E7£_ a - M£/j n) + Ffo, i n ) (15) 
.M being the difference operator defined by 



(16) 



m— 



Here and in the rest of the paper we adopt the convention that the summation 
is zero when the lower bound is larger than the upper bound. 
Equation (|15[) can be written in vectorial form as 



AU (n) =MU {n) +F {n) . 



(17) 
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Comparing this equation with SU = F one sees that 5 = A — M. The present 
method is implicit because one has to solve the tridiagonal system (fT5)l in order 
to get the numerical solution U n . Fortunately, the system can be efficiently 
solved by means of the Thomas algorithm since A is a strictly diagonally dom- 
inant matrix. For the case of constant timesteps, this method reduces to the 
one considered in [30j . 

The operator M. is a kind of difference operator with memory (which comes 
from the memory, the non-local character, of the fractional derivative) in the 
sense that its effect on U at time t n , U^ n \ depends on all the previous values 
E/W, • • • J/( n-1 '} = U{ n ~ x } . Where convenient, we will emphasize this 
fact by writing AiU^ as M. \U^- n ~ 1 ^ 1 \ as this last form has the virtue of making 
clear that the value of MXJ at time t n is obtained from the values of U at all 
the previous times t ra _i, • • • ,io- Then, the solution of (|17[) can be written as 

U {n) = A~ l M \u {n - x A + A^F^l (18) 

It is interesting to note that, if one discrctizes the Caputo derivative at time 
t n+ \ and the Laplacian at time t n , then one straightforwardly gets an explicit 
finite difference scheme 

n 

E f ^l + i [U? +l ~ U?] ~ S n [U? +1 2U? + U^} = F( Xj ,t n ). (19) 

m=0 

We will not explore this explicit method here because it is unstable if, for a 
given set of parameters 7, Ax, and K, the timesteps are not small enough. 
This indeed greatly reduces any advantage of a method with variable timesteps. 
Nonetheless, to examine how the non-uniformity of the timesteps affects the 
regions of stability in the space of parameters would be, from a theoretical 
perspective, an interesting topic in itself. For the case of constant timesteps, 



this explicit method becomes the one considered by the authors in [26 1 



3. Stability 

In order to check whether the finite difference scheme dU = F of equations 
(O or (TT5j) is stable, one studies how a perturbed solution U evolves with respect 
to the reference solution U, or, in other words, how (the size of) the perturbation 
v = U — U evolves in time. Usually, the initial perturbation is considered 
as the difference between the exact initial condition and its computer finite 
precision representation (round-off error). Because S is a linear operator, one 
sees that the perturbation satisfies Sv = 0, i.e., the same equation as U in ©, 
(USD, or (HP when F = 0: 

Av {n) =W B l (20) 

In order to analyze the stability of the present fractional diffusion difference 
algorithm, we use the von Neumann-Fourier technique introduced for this kind 
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of schemes in Refs. [21, 22, 23, 2^| : first one assumes that is described by 
a discrete Fourier series 

(21) 



A^i s 9 

9 



where the summation is carried over all the wave numbers q supported by the 
lattice, and then one analyzes the stability of the complete solution by analyzing 



The rationale for this 

(n) 



the stability of a generic Fourier q-mode, say £, 

procedure is that if any mode is stable then the complete solution vj l> written 
as a superposition of modes is stable too. 

Inserting the expression for the generic mode ^ n) e iqjAx into flU and using 
the definitions of A and M.U [see ([TBI) ] one S e ^ s 



(m) 
9 



m=0 



(22) 



where 



4 sin 



qAx 



(23) 



(7) _ 



1 and have defined f_i n = 



Here we have taken into account that 

Now we depart from our previous [23l . |24| quick method of analysis in terms 
of the so-called amplification factor (as it leads to sums involving Tm}n we have 
not been able to evaluate), and we use a procedure similar to the one followed 
by Murio in Ref. (3l|. However, it should be noted that our demonstration is 
more general (and difficult) as we do not assume, as in Ref. (3l|, that the (in 
general) complex quantity £^ is real and positive. 

and taking into account that Tm)n — T^_ t > for any temporal 
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mesh [see Appendix, equation (|A.8[) ]. one finds that 



n-1 

(1 + S n ) < - ^l,n) |4 m) | < |4 n_1 



ro=0 



n-1 

E 

m=0 



f 'f (t) _ 

\ m,n m— l,n 



=.(7) 
m — ! 

(24) 



where wc have defined 



t{n-l} 
So 



max • 



c(0) 

^9 



^9 



But srio 



(7) 



(7) 



1 [recall that „ = 1 and T 



t(n-l) 

^9 



■»9 — ^9 



; {n-l} 



(25) 
0] so that 
(26) 



because (1 + S n ) > 1. Obviously, J26J implies, a fortiori, that 



(„) 

9 



< 



(n-1) 
9 



(27) 
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for all n. Therefore, the perturbation of the generic mode remains smaller than 



or equal to its initial perturbation 



An) 



< 



AO) 



This means that the present 



difference method is unconditionally stable. 

Often the stability condition is expressed in terms of the 2-norm (or Eu- 



clidean norm) of the perturbation: |M n )|L = <-^x J^ - 



,(«) 



One can trans- 



form the stability (|2T[) to this form by means of Parseval's relation, 
2 , . 



,(»)! 



?9 



(») 



[46| . Using this expression (|27|) becomes 





< 


„(»-!) 




2 





(28) 



which implies that the norm of the perturbation remains bounded by the initial 
perturbation: ||tA n ) || < ||w < -°- ) || 2 . 

The solution of ([50)1 can be written as 



; (n) = A _1 M 



,{»-!} 



(29) 



so that the stability condition implied by (f2"B"]) means that the operator A 1 A4 
has the key property 



{„-!} 



< 



,("-!) 



(30) 



4. Consistency and stability implies convergence 

The basic property a difference scheme SU — F should have is that its 
approximate solutions should converge towards the exact solution of dU = F 
when the size of the spatiotemporal discretization goes to zero. In this case 
it is said that the method is convergent (46, 47 1 . In this section we show that 
for a whole class of fractional difference algorithms, similarly to the case for 
the usual difference algorithms, the consistency and stability of the difference 
scheme implies its convergence. (This result can be seen as the extension of 
the sufficient condition of Lax's equivalence theorem to fractional equations of 
the form du = F, or as the fractional counterpart of what is sometimes called 
fundamental theorem of finite difference methods |47j|.) 

(k) 

Let us define e) as the difference between the numerical and exact solution 



at the point (xj,t m ): e 



(k) 



U 



(fc) „,(fc) 



By definition du = 5u + R, SU = and 



du = 0, so that, because these operators are linear, Se = R. But 5 = A — Ai 
[see following equation (|T7|) ]. and then 



Ae {n) = M 



R (n) ^ e (n) = A -l M 



,{«-!} 



A- 1 ^. (31) 



Therefore 



< 



A~ l M 



= {«-!} 



(32) 
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and, because the method is stable, (|30p holds and (|32p becomes 



e <») 


< 


e ("-D 


4 


Ik" 1 !!, 






2 




2 





(33) 



Then, by induction, one finds that 



< A' 



E ||^ (Tn) 



m.—l 



(34) 



If, when the size of the discretization mesh goes to zero, (i) ||-R^ Tn ' 1 || 2 S oes to 
zero (i.e., if the method is consistent) and (ii) ||^4 _1 || 2 is bounded, then (|34[) 
implies that the error ||e^ n '|| 2 goes to zero, i.e., the method is convergent. In 
particular, the finite difference method of this paper is convergent because, (i) 
the method is consistent [see following equation ©] and (ii) ||-4 _1 || 2 is always 
bounded since A is a real symmetric tridiagonal Toeplitz matrix with all its 
eigenvalues greater than 1 [47[. Note that the present demonstration is quite 
general as only generic properties of the operators Ai and A appearing in the 
finite difference equation (linearity, boundedness, stability) are required. In fact 
this demonstration is parallel to the standard procedure to prove this result for 
non- fractional diffusion problems [48| . This way, one realizes that some aspects 
of fractional difference methods can be included into the realm of standard non- 
fractional finite difference theory. In fact, fractional difference methods could 
be seen as particular cases of multilevel schemes where the number of levels 
increases with time. 



5. Dispersion of a flux of subdiffusive particles: Numerical results 

In this section we test the finite difference scheme (|9|) by solving a problem 
that describes the dispersion of a constant flux of anomalous subdiffusive par- 
ticles appearing at a given place. In particular, in our problem the flux is 1 as 
one particle is released every unit time at times t = 0,1,2... from a source of 
particles placed at x — in a one-dimensional infinite medium. We shall assume 
that the probability density u(x, t) of finding the particle at position x at time t 
follows the fractional diffusion equation du = 0. The release of a particle at time 
n is described by the introduction of the Dirac delta function 8{x) at that time 
[see following equation <j2j) ] - The probability distribution associated with each 
particle is just the propagator (or Green's function) G(x, t — k) of the problem, 
i.e., the solution of du = in the unbounded space when there is only a single 
Dirac delta function that appears at t = k. This propagator can be explicitly 
written in terms of Fox's H function 0, |49| : 



G(x,t) 



rrlO 
H ll 



(1-7/2,7/2) 
(0,1) 



(35) 
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Finally, given the linear nature of the problem, its exact solution is just a su- 
perposition of Green functions: 

L*J 

u(x,t) =^TG(x,t-k) (36) 

fe=0 

where [^J is the floor function. In the formalism of ((1), the function F of this 
problem is 

LtJ 

F(x,t) =5(x)J2$ h Ht-k) (37) 

k=0 

with 

*M( t -fc) = i!le(f-fc), (38) 

9 being the Heaviside step function. 

The problem just presented, with the added feature that the particles may 
react and disappear, describes the development of morphogen gradients (a key 
concept in developmental biology) when the anomalous diffusion of the mor- 
phogens is accounted for by means of a CTRW model [Ho| • Because knowledge 
of the long-time distribution of morphogcns in the medium is especially impor- 
tant, numerical methods that can use large timesteps (to go fast in time) and 
small timesteps (to take account of regions of large variability when the particles 
are released) would seem to be particularly convenient. 

Regarding the boundary conditions, although it is of course not possible 
to use boundary conditions at infinity in numerical calculations, one can solve 
the problem by means of some convenient boundary conditions, say absorbing 
boundary conditions u(x = —L/2,t) = u(x = L/2,t) = 0, so far away from the 
source at x = that, for a finite t, their impact on the solution is negligible in 
the region of interest which is, by construction, far away from the boundaries. 
In our computations, where we have chosen K = 1 and 7=1/2 and the largest 
value of time is t = 2, we used L — 20. For this choice the exact value of 
the solution at the boundaries is u(±10, 2) ~ 5 x 10~ 5 , which can be safely 
approximated by in our numerical calculations (as we will show shortly, the 
errors are well above 10 -4 ). The nodes of the spatial mesh we used were placed 
at x 3 = j Ax with j = -N/2, -N/2 + 1, ■ • • , N/2, N = 100, and Ax- = 0.2. 
Finally, every Dirac delta function 6(x) was approximated by its discrctized 
version: 

x, \ / j = 0, , . 

S W-\0, j?0. (39) 

Note that the discontinuity of the solution due to the periodic appearance of 
Dirac delta functions at times t = k (k integer) has yet to be worked out. Al- 
though in section[2]we only consider the case in which the solution is continuous, 
discontinuities can be easily handled by our algorithm taking into account ([7|. 
which is valid even for discontinuous functions at the discretization times t m . 
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This way, keeping the symbol Uj for the estimate of u(xj,t m ) and denoting 
by Vj the estimate of u{xj,tf n ), ([9]) becomes 

n-l 

E S [ f/ r +1 - v n ~ w+i - 2c/ ; + ^--i] - ^> **)• ( 4 °) 

m=0 

The algorithm in this case can be written as in (|17j) but now with 

n-2 

MU { / 1) = VJ 1 - 1 - E ^iU [^ m+1 - V™] • ( 41 ) 

m— 

In particular, the appearance of a Dirac delta function at x = at time t = k 
is handled by simply making Vq = 1/ Ax + Uq whenever t n = k. 

In Fig. Q]we compare exact solutions and numerical solutions obtained by 
means of the present implicit method when the non-uniform discretization is 
given by this adaptive timestep: 

Wi - t m = min [1(T 4 coth [\g(t m )\/1000] , 0.02] . (42) 

Here g(t m ) = (UH\ - 2£/g" + U? l )/(Ax) 2 is the three-point centred estimate of 
the spatial second derivative at x = and time t m : g{t m ) ~ d 2 u(x, t m )/dx 2 \ x =o, 
i.e., an estimate of the curvature of the solution at the place where it changes 
most abruptly. The agreement is excellent. In fact, even at x = 0, the place 
where the numerical solution is worst, the results are quite good (see Fig. [5]). 
The relatively large error for short times after the appearance of the Dirac delta 
function is due to its crude approximation by 

Although formula (pf2"|) for the temporal spacing t m+ i~t m is indeed arbitrary, 
it has some convenient features: it has prefixed minimum and maximum values 
(0.0001 and 0.02, respectively) and t m+ i — t m increases (decreases) when the 
curvature at x = 0, d 2 u(x 1 t m ) / dx 2 \ x= $ decreases (increases). This guarantees a 
more thorough description of the solution where necessary, namely, those regions 
where the solution changes most quickly. For example, when (|42l) is used, 5 
timesteps are employed in the tiny temporal region between t = and t = 0.01 
(a region where the solution changes abruptly; see Fig. [2]); however, only 59 more 
timesteps are necessary to cover the comparatively much larger region between 
t = 0.01 and t = 1. Note that this makes the procedure with variable timesteps 
as least as accurate as the method with fixed timesteps i m +i — t m — 0.001 in 
the difficult region near t = (sec Fig. [3]). However this is not reached at at 
the expense of computational cost: compare the 64 steps [17 steps] required 
by the variable timesteps method to reach t = 1 [t = 0.1] with the 1000 steps 
[100 steps] required when £ TO +i — t m = 0.001. This difference is far from minor 
because, from ([5]), one sees that the computation time required to calculate 
the solution after n timesteps increases roughly as n 2 . In fact, in our actual 
calculations we found that the CPU time spent to evaluate the solution at t = 1 
with t m+ i — t m = 0.001 was well above one hundred times! longer than when 
(|4"2")l was used. 
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X 



Figure 1: Exact solution u(x, t) (lines) and numerical solutions Uj with variable timesteps 
given by ( 1421 ) (symbols) for the problem described in the main text with (from top to bottom 
at x = 0) t = 1.0004 (66 timesteps, down triangles), t = 4.08 X 10 — 4 (one timestep, squares), 
t = 0.034 (10 timesteps, diamonds), t = 2.0 (141 timesteps, circles), and t = 1 (65 timesteps, 
up triangles). The lines for t > 1 are dashed. Inset: logarithmic scale and tails of the solutions. 
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Figure 2: Exact solution u(0, t) (line) and numerical solution Uq 1 ^ (symbols) at the origin 
x — for the problem described in the main text for fixed timesteps with i m +i — t m = 0.001 
(circles) and for variable timesteps given by 1142 I I (solid squares) in the time interval < t < 0.1. 
Inset: solutions in the whole time interval. In this panel only 1 of every 20 points for the case 
with fixed timesteps is plotted. 
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Figure 3: Absolute difference between the exact solution and the numerical solution at x = 0, 
e(0, t)|, for the problem described in the main text for fixed timesteps with t m ^i —t m = 0.001 
(circles) and variable timesteps given by (142 [ i (squares) in the time interval < f < 0.1. Inset: 
error in the whole time interval. In this small panel only 1 of every 20 points for the case with 
fixed timesteps is plotted. 



6. Conclusions and final remarks 

An adaptive finite difference method with non-uniform timesteps for frac- 
tional diffusion equations has been described. The method allows the temporal 
mesh to be adapted to the behaviour of the solution by sampling thoroughly 
those regions where the variation of the solution is large -so as to maintain 
the accuracy of the method- but sparsely wherever the solution varies mildly 
-and so advancing fast in time and reducing the number of timesteps needed 
in the computation. This last feature is especially convenient when numerically 
solving fractional problems since the computational cost in this case scales as 
the square of the number of timesteps. Compare this with the typical linear 
dependence for non-fractional diffusion problems. 

The finite difference method discussed in this paper was obtained from the 
fractional diffusion equation in the Caputo form by discretizing (i) the Lapla- 
cian by means of the standard three-point centred formula and (ii) the Caputo 
derivative with a generalization of the LI formula to non-uniform meshes. For 
the case of uniform timesteps, the method becomes the one considered by Liu 



et al. [30| and Murio [31|. The corresponding truncation error is of second 
order in the size of the spatial mesh and of first order in the size of the tem- 
poral mesh. A von-Neumann stability analysis was employed to prove that the 
method is unconditionally stable. Also it was shown in a quite general form 
that the consistency and stability of the method implies its convergence. In 
fact, this implication, sometimes called the fundamental theorem of finite dif- 
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ference methods, it is not restricted to the particular finite difference scheme 
considered in this paper, but it is a quite general result that depends only on 
some general properties of the fractional finite difference scheme. 

The present method can be generalized and extended in several ways. First, 
one could employ other different non-uniform discretizations of the fractional 
derivative, say of higher-order accuracy, or even consider non-uniform discretiza- 
tion in the space. Also one could explore criteria or formulas for efficiently choos- 
ing the size of the timesteps. Finally, another interesting avenue to explore is 
the application of adaptive methods to other kinds of fractional equation, such 
as fractional diffusion-wave equations, fractional Fokkcr-Planck equations, and 
fractional equations with Riemann-Liouville derivatives. 
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Appendix A. Non-uniform LI formula for the Caputo derivative 

We start by rewriting the Caputo derivative defined in @: 



dr> y (t) 



dt~i 



1 



r(i- 7 ) 



n-l r t r 

E 



- 1 dy(r) 



(t - rp dr 



(A.1) 



*-*n v ' ' m =0 * 

with t = and n > 1. Proceeding as is standard in the so-called LI method [l[ 
the ordinary derivative inside the integral with lower limit t m is approximated 

by 

dy(t) _ y(t m+1 ) - y(t+) 



dt 



+ 0(t m +x — £ m ), t m < t < t 



771+ 1 



(A.2) 



so that 



dV 



1 „ 



E 



y(t~ +1 ) - y(t+ ) 



t-t r(i — 7) t m+ i - 1,- 



dr(t - t) _t + R(t n ) 

(A.3) 



where y(t m ) and y(t^) are the values of y(t) when t goes to t m from the left 
and right, respectively, and 



R(t n ) 



m=0 



dr{t - t)~ 7 



(A.4) 



is the temporal truncation error. The evaluation of the integral of (|A.3|) is 
straightforward and this equation becomes 



dty(t) 



dt~> 



1 n — 1 

E T ™U - !/(*£)] + Rt{n) (A.5) 



r(2 -7) 
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where T^} n is given is (|13[) . Therefore, from (|A.4[) one sees that 

\R(t n )\ < rn 1 > mgx |i? m | / " dT{t-T)-i (A.6) 

1 (^1 — 7j 0<m<rt — 1 Jq 

so that Rt(n) < Cmaxo< m <„-i |i? m |, C being a constant of order i 1-7 , i.e., 



max (i m+ i - ^ 

0<m<n-l 



(A.7) 



Finally, we shall prove that 

t(t) -> t^t) f a Si 

-^771,71 ^ -^m— l,n' ^xi.uy 

a result we used in (j2"4")l . From the definition of T^„ and making the change 
z m = tn — t m , one sees that 

T<& = 4r7 _~ Z "' +71 . (A.9) 

Zm 

Let us define f(z) = z 1 ^ 7 . To prove (|A.8[) is equivalent to proving that 

/(Zm+l) - /(Zm) < f(z m ) - f(z m -l) ^ 
Z'tn+1 Zm Z m Z m _i 

which is true for any monotonous growing function f(z) with monotonous de- 
creasing derivative f'{z) such as f(z) = z 1-7 . 
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